#Code is written for R Open 3.5.1, which can be downloaded here: https://mran.microsoft.com/release-history
#The CRAN mirror snapshot was taken on 2018-08-01.  R Open 3.5.1 automatically installs packages from this snapshot.

library(data.table) #setnames function, fread

options(max.print = 2000)

# load ONE of these two .csv's.
    # In the first, the "search" variables are derived with six terms: smog + no2 + air pollution + ozone + pm2.5 + AQI
    # In the second, the "search" variables are derived with five terms: smog + no2 + air pollution + ozone + pm2.5
    # The Table 1 results shown in the manuscript use the us_pollution_paper_cbsas_dfmain (google 6 terms).csv data.

df.main <- fread("./us_pollution_paper_cbsas_dfmain (google 6 terms).csv")
#df.main <- fread("./us_pollution_paper_cbsas_dfmain (google 5 terms).csv")

# Public concern / Google Search Index models
    m1 <- lm(search.adj ~
        search.adj.tm1 +
        no2.change +
        all_enf.change +
        no2.state.change +
        all_enf.state.change +
        vote.prop +
        #factor(year) +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))      
    m2 <- lm(search.adj ~
        search.adj.tm1 +
        no2.state.change +
        all_enf.state.change +
        vote.prop +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m3 <- lm(search.adj ~
        search.adj.tm1 +
        no2.change +
        all_enf.change +
        vote.prop + 
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))

    summary(m1)
    summary(m2)
    summary(m3)
        

# Policy / CAA enforcement action models
    m1 <- lm(all_enf ~
        all_enf.tm1 +
        search.adj.change.tm1 +
        no2.change.tm1 +
        search.state.adj.change.tm1 +
        no2.state.change.tm1 +
        #factor(year) +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m2 <- lm(all_enf ~
        all_enf.tm1 +
        search.state.adj.change.tm1 +
        no2.state.change.tm1 +
        #factor(year) +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m3 <- lm(all_enf ~
        all_enf.tm1 +
        search.adj.change.tm1 +
        no2.change.tm1 +
        #factor(year) +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))

    summary(m1)
    summary(m2)
    summary(m3)

# Outcomes / NO2 models
    m1 <- lm(no2 ~
        no2.tm1 +
        grp.change +
        search.adj.change.tm1 +
        all_enf.change.tm1 +
        search.state.adj.change.tm1 +
        all_enf.state.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m2 <- lm(no2 ~
        no2.tm1 +
        grp.change +
        search.state.adj.change.tm1 +
        all_enf.state.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m3 <- lm(no2 ~
        no2.tm1 +
        grp.change +
        search.adj.change.tm1 +
        all_enf.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))

    summary(m1)
    summary(m2)
    summary(m3)

# Outcomes / PM2.5 models
    m1 <- lm(pm25 ~
        pm25.tm1 +
        grp.change +
        search.adj.change.tm1 +
        all_enf.change.tm1 +
        search.state.adj.change.tm1 +
        all_enf.state.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
        summary(m1)
    m2 <- lm(pm25 ~
        pm25.tm1 +
        grp.change +
        search.state.adj.change.tm1 +
        all_enf.state.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))
    m3 <- lm(pm25 ~
        pm25.tm1 +
        grp.change +
        search.adj.change.tm1 +
        all_enf.change.tm1 +
        name.state
        , data = df.main, subset = (year > 2009 & memi == 1))

    summary(m1)
    summary(m2)
    summary(m3)
